source("/home/jc227089/evo-dispersal/KBGrad/KBGradfunctions.R")
setwd("/home/jc227089/Temp/")

stable.domain<-50

K<-30
lmda<-10
b<-0
n<-2*stable.domain*K
spX<-stable.domain
ngens<-100
Hmean<-0
Dmean<-0
h2H<-0.5
h2D<-h2H
VPH<-0.25
VPD<-0.05


##############################
mother<-function(n, spX, K, lambda, b, ngens, Hmean, Dmean, h2H, h2D, VPH, VPD, evovar, stable.domain=spX){
  pop<-init.inds(n=n, spX=spX, Hmean=Hmean, Dmean=Dmean,
	h2H=h2H, h2D=h2D, VPH=VPH, VPD=VPD)
	ztrace<-c()
  for(ii in 1:ngens){
    pop<-repro.disp(popmatrix=pop, K=K, lambda=lambda, b=b, 
    	h2H=h2H, h2D=h2D, VPH=VPH, VPD=VPD, evovar=evovar, stable.domain=stable.domain)
    save(pop, file=paste("out", ii, ".RData", sep=""))
  }
}

mother(n=n, spX=spX, K=K, lambda=lmda, b=b, ngens=ngens, Hmean=Hmean, 
		Dmean=Dmean, h2H=h2H, h2D=h2D, VPH=VPH, VPD=VPH, evovar=TRUE, 
		stable.domain=stable.domain)
		
##############################


# plots population density and mean trait values through space
plotter.mean<-function(popmatrix, filename=NULL, Nx.r, H.r, D.r, ...){
	popmatrix<-popmatrix[order(popmatrix[,"X"]),]
	summ<-sum.metrics(popmatrix, TRUE)
  	if (!is.null(filename)) pdf(file=filename, width=5, height=10)
	par(mfrow=c(3, 1), mar=c(5,5,2,2), cex.lab=1.4)
	
	plot(summ[,"X"], summ[,"Nx"], xlab="Location", ylab="Population density", type="l", ylim=Nx.r, bty="l", ...)
	
  	plot(summ[,"X"], summ[,"MeanH"], xlab="Location", ylab="Mean trait value", pch=19, col="darkorange", 
  		type="l", ylim=H.r, bty="l", ...)
	
	plot(summ[,"X"], summ[,"MeanD"], xlab="Location", ylab="Mean dispersal value", pch=19, col="darkgreen",
		 type="l", ylim=D.r, bty="l", ...)
	
	if (!is.null(filename)) dev.off()
}


poplist<-vector(length=ngens, mode='list')
max.x<-0
min.x<-0
min.Nx<-0
max.Nx<-22
min.D<-(-1)
max.D<-2
min.H<-0
max.H<-0.5
for (ii in 1:ngens){
	load(file=paste("out", ii, ".RData", sep=""))
	poplist[[ii]]<-pop
	max.x<-max(max.x, max(pop[,"X"]))
	min.x<-min(min.x, min(pop[,"X"]))
	#max.D<-max(max.D, max(pop[,"D"]))
	#min.D<-min(min.D, min(pop[,"D"]))
	#max.H<-max(max.H, max(pop[,"H"]))
	#min.H<-min(min.H, min(pop[,"H"]))
}

pdf("expansion.pdf", width=5, height=10)
for (ii in 1:ngens){
	plotter.mean(popmatrix=poplist[[ii]], xlim=c(min.x, max.x), Nx.r=c(min.Nx, max.Nx),
		H.r=c(min.H, max.H), D.r=c(min.D, max.D))
}
dev.off()